Brittle Crack Roughness in Three-Dimensional Beam Lattices 
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The roughness exponent is reported in numerical simulations with a three dimensional elastic 
beam lattice. Two different types of disorder have been used to generate the breaking thresholds, 
i.e., distributions with a tail towards either strong or weak beams. Beyond the weak disorder regime 
a universal exponent of £ = 0.59 ± 0.01 is obtained. This is within the range 0.4 < £ < 0.6 
reported experimentally for small scale quasi-static fracture, as would be expected for media with 
a characteristic length scale. 
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The interest in crack morphology can be traced back 
to Mandelbrot et al., who introduced a mathematical 
framework for describing rough surfaces in terms of frac- 
tal geometry More specifically, crack surfaces have 
been shown to be self-affine in the sense that they sat- 
isfy certain scaling laws, where a characteristic exponent 
governs the asymptotic behaviour [2j • Since the quantita- 
tive properties of crack surfaces have been made readily 
accessible through the development of modern imaging 
techniques, predictions of computer models [3j can be 
compared with the results of practical experiments to 
verify or discard theoretical assumptions. 

Much of the theoretical work done so far has been 
based on the random fuse model i.e., a regular lattice 
of conducting elements which irreversibly burn out once 
their individual thresholds have been exceeded. Break- 
down is driven by a potential difference between two op- 
posing boundaries and the analogy of Kirchhoff 's equa- 
tions with linear elasticity is the reason why this model is 
referred to as a scalar model of fracture. The result ob- 
tained for the roughness exponent in simulations with the 
fuse model is C = 0.74(2) in two dimensions 0- Recently, 
a somewhat different result, £ — 0.84(3), has been re- 
ported . Experimental studies in two dimensions have 
been carried out by Poirer et al. @, who considered a 
stacking of collapsible cylinders (drinking straws) to ob- 
tain £ = 0.73(7), by Kertesz et al. |8j for tear lines in 
wet paper, obtaining ( » 0.73, and by Eng0y et al. 
for crack lines in thin wood plates, yielding a roughness 
exponent of £ = 0.68(4). Results also differ according 
to the dimensionality, i.e., in three dimensions the fuse 
model result is C = 0.62(5) [HJ]. This result, however, 
is very different from the universal value suggested by 
Bouchaud et al. 0, i.e., £ = 0.8, a value which has been 
confirmed in experimental studies Hence the appar- 
ent agreement of the fuse model with experiment in two 
dimensions is a coincidence. 

An important point to be considered with the fuse 
model, however, is the fact that it models electrical 



breakdown rather than brittle elastic fracture. A dif- 
ferent model which incorporates elasticity is the Born- 
model yj|]. Here the material is modeled as a network 
of elastic springs, each spring being free to rotate at its 
ends. In two dimensions £ = 0.64(5) is obtained with 
this model [lif. As with the fuse model, the exponent 
drops when going from two to three dimensions. The re- 
sult, £ = 0.5 [l5j, is also significantly different from the 
universal value of £ = 0.8. Recent findings, however, in- 
dicate that there should exist two different regimes for 
the scaling behaviour of cracks. As was first shown by 
Milman et al. [lfij . a much lower exponent £ « 0.5 is 
found on small length scales. Daguier et al. identified 
two self-affine fracture regimes, where £ ~ 0.5 is associ- 
ated with small length scales and slow crack growth and 
£ = 0.8 involves dynamic effects and is associated with 
larger length scales . 

A problem with the Born model is that it is not rota- 
tionally invariant 0] . A different model which does not 
suffer from this drawback is the elastic beam model • 
Here the beams are rigidly connected at the nodes so 
as to preserve the angle between any two neighbouring 
beams. Rotations thus induce flexing and twisting de- 
formations, while linear displacements induce transverse 
shear and axial forces. The exponent obtained with the 
beam model is £ = 0.86(3) in two dimensions j2lj, an d 
shows that the scalar and vectorial descriptions of frac- 
ture need not necessarily produce similar results. Since 
the beam model provides a more realistic description of 
brittle elastic fracture than does either the fuse model 
or the Born model it is of great interest to see what the 
result is in three dimensions. 

The lattice used in our calculations is a regular cube of 
size L x L x L, where each node is connected to its near- 
est neighbours by linearly elastic beams of unit length. 
Forces acting on the nodes have been deduced from the 
effect a concentrated end-load has on a beam with no 
end-restraints [21j . A coordinate system is placed on each 
node, and the enumeration of the connecting beams fol- 
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FIG. 1: Enumeration scheme on the beams of a cube lattice 
connecting node i to its j = 1-6 nearest neighbours, showing 
the coordinate system with i as its origin. 

lows an anti-clockwise scheme within the XY-plane, i.e., 
beginning with the beam which lies along the positive X- 
axis and ending with that which extends upwards along 
the positive Z-axis, see Fig. ^ 

At each stage of the breaking process, the updated 
displacements for each node is obtained from 
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which is solved iteratively via relaxation using the conju- 
gate gradient method. This minimizes the elastic energy 
to obtain those displacements for which the sum of forces 
and moments on each node vanish, i.e. the mechanical 
equilibrium. 

Six terms contribute to each of the force components 
in Eq. (JJJ. With the neighbouring nodes fixed, a transla- 
tional displacement 5x = Xj — Xi of node i induces axial 
forces in beams 1 and 3, and transverse forces in beams 



2, 4, 5 and 6. With the expressions for force and moment 
defined in Refs. [i(J and j^, this gives 

X i = x Af ) + x Af ) + A j \ (2) 

with similar expressions for Yi and Z{. Likewise, a ro- 
tational displacement Su = Uj — Ui about the X-axis 
causes a torsional moment in beams 1 and 3, and flexure 
in beams 2, 4, 5 and 6. Hence, 

Ui = x TP+ x Tf+ £ X M®, (3) 

now with similar expressions for Vi and Wi. 

To express the thirty-six force components in Eq. 
more compactly, 

^■ = nW (4) 

n=0 

and 

Sj = ( D'o (5) 

arc defined for notational convenience, to keep track of 
the signs. The Kronecker delta, moreover, has been used 
to construct 

X<f t l=5 sj + 6 tj (6) 

and 

7&}=5 kj (l-6.j)(l-8 tj ), (7) 

operators which include or exclude s and t from the sum 
over neighbouring beams. 

Eq. (J2J can then be stated on the form 



(8) 



where, exchanging indices (1,3) and (2,4), and interchanging v with u and x with y, Yi is obtained from Xj. In the 
same way, Zi is obtained from Xi by an exchange of the indices (1,3) and (5,6), by interchanging w with u and x with 
z, and by replacing rj with —Sj. 

Next, Eq. © can be stated on a similar form, i.e., 
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where Vi, for angular displacements about the K-axis, angular displacements about the Z-axis, follows from Ui 
is obtained from Ui by the same substitutions as in by use of the Xi — > Zi substitutions, now however, with- 
Xi — > Yj. Likewise, the expression for Wj, relevant to out changing rj to — Sj. Instead, the change from z to x 



3 



is made so that — Sx replaces Sz. 

Breaking thresholds are generating by raising a ran- 
dom number r on the interval [0, 1] to a power D. The 
magnitude of D then corresponds to the strength of dis- 
order, with the distribution having a tail towards weak 
beams for D > and towards strong beams for D < 0. 
This is the most simple way of incorporating scale invari- 
ance, i.e., those mechanisms responsible for the multifrac- 
tal behaviour of disordered systems, into the distribution 
of thresholds [2^]. For each sample two casts are gener- 
ated for the thresholds, one for the flexural strength and 
one for the axial strength. 

Fracture is initiated by imposing a uniform displace- 
ment on all nodes defining the top surface of the cube. 
The first beam to break is the one with the weakest axial 
strength, whereupon the location of subsequent breaks 
depends on a complex interplay between quenched disor- 
der and the constantly evolving non-uniform stress-field. 
Each time a beam is removed from the lattice, Eq. (JJJ is 
used to obtain the new equilibrium displacements. This 
is continued until a surface appears which divides the 
sample in two separate parts, see Fig. |5J 

Results obtained for D = 1, 2 and 4 are shown in 
Figs. 01 El and respectively, where the roughness W has 
been calculated for systems ranging in size from L = 3 
to L = 40. Disregarding finite-size effects for small L, 
the relationship between L and W is clearly self-affine. 
Hence, defining the roughness as the root-mean-square 
of the variance perpendicular to the (average) fracture 




FIG. 2: Fracture surfaces, showing the lower remaining part 
of a cube lattice of size L = 40 after it has broken completely. 
Four samples are included for D = 1 on the left, for D = 2 at 
center, and for D = —4 on the right. 
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FIG. 3: Logarithmic plot, showing the roughness, W , as a 
function of the system size, L, for disorder D = 1. Filled 
circles denote those points to which the straight line, with 
slope £ = 0.77(1), has been fit. 
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we have W ~ L^, with the results of W x and W y being 
statistically the same. 

The exponent in Fig. is C = 0.77(1), and does not 
belong to the self-affine regime of fracture. This is be- 
cause D = 1 in a three dimensional lattice generates too 
little disorder. The qualitative difference between the 
fracture surfaces obtained for D — 1 and higher values 
of D is shown in Fig. [3 A universal exponent emerges 
as the disorder is increased, however, with the exponents 
obtained for D = 2 and D = 4 being C = 0.62(2) and 
£ = 0.592(5), respectively. 

The result also seems to hold for disorders with a tail 
towards strong beams. For D = —4 the exponent is 
C = 0.65(2), see Fig. El Although this is slightly different 
from the best value obtained for D > 0, which we take 
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FIG. 4: Logarithmic plot, showing the roughness, IF, as a 
function of the system size, L, for disorder D = 2. Filled 
circles denote those points to which the straight line, with 
slope £ = 0.62(2), has been fit. 
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FIG. 5: Logarithmic plot, showing the roughness, W , as a 
function of the system size, L, for disorder D = 4. Filled 
circles denote those points to which the straight line, with 
slope £ = 0.59(1), has been fit. 

to be C = 0.59(1), the small discrepancy is probably due 
to the fact that D = —4 lies just within the transient 
regime towards strong disorders for distributions of type 
D < 0. In two dimensions the D < stress-strain curve 
displays a cross-over towards stable crack-growth when 
\D\ ~ 2. This cross-over value is probably higher in the 
cube lattice, due to the stronger constraints imposed on 
the crack path in three dimensions. 

To summarize, the roughness exponent £ = 0.59(1) is 
obtained numerically for a cubic lattice of elastic beams. 
This is lower than the experimental value relevant to 
large scales, £ = 0.8, but within the range reported for 
small scales, 0.4 < ( < 0.6. This lends further evi- 
dence to the existence of a different scaling regime for 
slow crack growth involving characteristic, or "small", 
length scales |2^| . The agreement of our result with the 
exponents reported in this regime stems from the dis- 
cretization in terms of a beam lattice, i.e., a system which 
involves a characteristic intermediate length scale, as in 
the case of Cosserat elasticity. 
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FIG. 6: Logarithmic plot, showing the roughness, W , as a 
function of the system size, L, for disorder D = —4. Filled 
circles denote those points to which the straight line, with 
slope ( = 0.65(2), has been fit. 
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